Dynamic response functions for the Holstein-Hubbard model 
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We present results on the dynamical correlation functions of the particle-hole symmetric Holstein- 
Hubbard model at zero temperature, calculated using the dynamical mean field theory which is 
solved by the numerical renormalization group method. We clarify the competing influences of the 
electron-electron and electron-phonon interactions particularity at the different metal to insulator 
transitions. The Coulomb repulsion is found to dominate the behaviour in large parts of the metallic 
regime. By suppressing charge fluctuations, it effectively decouples electrons from phonons. The 
phonon propagator shows a characteristic softening near the metal to bipolaronic transition but 
there is very little softening on the approach to the Mott transition. 

PACS numbers: 71.10.Fd,71.30.+h,71.38.-k 
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I. INTRODUCTION 

Phonons are important in metallic systems. They af- 
fect the electronic behaviour in diverse ways causing the 
dominant temperature dependent contribution to the re- 
sistivity, an enhancement of the specific heat, and an 
attractive effective electron-electron interaction, which 
may induce superconductivity. One can expect strong 
electron-phonon interactions in strongly correlated sys- 
tems, such as heavy fermion compounds 1 , where the ra- 
dius of the rare earth ions is very sensitive to the /- 
electron occupation. The coupling of electronic states to 
the lattice also plays an important role in the anomalous 
electronic behaviour of the manganites 2-4 . The super- 
conductivity in the fullerides 5 appears to be induced by 
electron-phonon interactions within a strongly correlated 
electron band. 

However, the fully quantum mechanical treatment of 
lattice effects in these compounds has so far received 
comparatively little attention. The reason for that lies 
in the difficulty in handling strong electron-phonon cou- 
pling together with strong electron-electron repulsion. 
Recent advances in non-perturbative techniques, most 
notably the development of the dynamical mean field the- 
ory (DMFT) 6 , enable one to investigate the interplay of 
electron and phonon interactions in these systems. 

One of the most prominent models in the field of 
strongly correlated electron systems is the Hubbard 
model 7 which been used extensively to study the effects 
of strong local electron-electron interactions. Especially 
within the framework of DMFT 6 the nature of the Mott 
transition could be clarified 8,9 . The Holstein model 10 
has been used to study polaronic effects in the absence 
of electron-electron repulsion, mainly in the limit of low 
electron density. Systems with finite electron density 
have received less attention, but have been studied within 
the DMFT using Monte Carlo 11 , exact diagonalization 
(ED) 12 , perturbative approaches 13-15 and the numerical 
renormalization group (NRG) 16 . 

Even less studied has been the interplay of 
electron-electron interactions and the electron-phonon 
coupling 17-21 . The most natural model incorporating 




FIG. 1: Zero temperature phase diagram of the half filled 
Holstein-Hubbard model for ujq = 0.2 and a semielliptic band 
of width W — 4 as calculated in Ref. 22. The dashed line 
is the locus of U c s = 0. Points on the dot-dashed line fulfil 
(n T n x ) = 1/4. 



both these terms is the Holstein-Hubbard model defined 
by the Hamiltonian 
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Here U describes the electron-electron interaction within 
a band of dispersion e(fc). The electron density + 
mi at site i couples linearly to the local displacement 
operator Xi = (b i + 6j)/^/2mojo with an electron-phonon 
coupling g. The phonons are assumed to be dispersionless 
(local Einstein phonons) with energy uj and m is the 
mass of the vibrating ions. 

Recently, the T = phase diagram of the half filled 
Holstein-Hubbard model has been calculated 22,23 using 
various local approximations, and is shown in Fig. 1. 
The results are for a semi-elliptic band of width W = 4 



2 



and phonon frequency ujq — 0.2 (m = 1). All types of 
long-range order are excluded. We can distinguish three 
different phases: metallic, bipolaronic and Mott insulat- 
ing phase. The metallic phase is always found to be a 
Fermi liquid. The latter two phases display a gap in the 
one-electron spectra and will be referred to as gapped 
phases. The parameters of Fig. 1, namely W — 4, m = 1 
and u) — 0.2 are taken throughout this paper, unless men- 
tioned otherwise. 

One way of simplifying the discussion of the Holstein- 
Hubbard model would be to integrate out the phonons 
in order to derive an effective model for the electrons. 
The effective action is then governed by the w-dependent 
potential 24 

, .. 2 a 2 LUn 

U eS (co) = U+-^°. (2 ) 

This reduces to the static quantity U e g = U — 2 g 2 /ujq 
in the limits u> — > and loq — > oo. Therefore, on small 
energy scales, we expect an effective attractive or repul- 
sive interaction depending on the value of U e ff. For small 
phonon energies luq, however, there is no reason to expect 
the static quantity J7 e ff to be sufficient to characterize the 
behaviour of the Holstcin- Hubbard model. The proper- 
ties of this model are expected to depend on both U and 
g individually. This is also apparent from the phase di- 
agram in Fig. 1. The dot-dashed line, where the double 
occupancy takes the value of a free system (n^ni) = 1/4 
is quite distinct from the line U e s = 0. The aim of this 
study is to clarify the interplay of these different types of 
interactions. 

An interesting question about this model is whether 
the suppression of local charge fluctuations by a signif- 
icant Hubbard U will effectively decouple electrons and 
phonons on all energy scales. It is possible, however, 
that for the behaviour on low energy scales the effective 
onsite-t/ is strongly renormalized and hence the electron 
phonon coupling will have a comparatively stronger effect 
on energy scales ui < ujo (see Ref 17,20). 

In this paper, we present results on dynamical elec- 
tronic and phonon response functions for the Holstein- 
Hubbard model as calculated within the dynamical mean 
field theory 6 (DMFT). The effective impurity problem is 
solved using Wilson's numerical rcnormalization group 
(NRG) method 25 extended to treat the electron-phonon 
coupling 24 . The combination of DMFT and NRG has 
proven to be a reliable and effective method for calcula- 
tions of both metallic and insulating phases at T = 0. It 
has been applied to the Mott transition in the pure Hub- 
bard model 26 and to the pure Holstein model 16 . In our 
calculations we use the discretization parameter A = 1 .8 
and retain up to 1200 states. The value of A = 1.8 gives 
rise to a bandwidth correction factor Aa = 1.029 (see Eq. 
(5.42) in 25) which we take into account in the NRG pro- 
cedure. For the calculation of dynamical response func- 
tions we use the method described in Ref. 27. 

Apart from the usual local one-electron Green's func- 
tion G a (uj) = ((c it7 ; c}o-)) w we wm b e interested in two 



different phonon Green's functions, 

d{u) = ((b i ;bt)) u , (3) 
D{u>) = 2 Wo (( Xi ; Xi )) u = {{h + b\ ■ b t + b\)) u . 

Since we will treat only homogeneous phases here and 
exclude all types of long-range order, we have dropped 
the site indices in all of the definitions. The spectral 
weight of d(uj) will be denoted by Pd{w) = —\ lmd(uj + 
i0 + ), and similar for all other Green's functions. The 
dynamical charge and spin susceptibility are given by 

X cM = «0;0» w , (4) 
X a (u) = «n iT - n n ;n iT - n il )) u , 

respectively. The operator O represents the electronic 
term that couples to the phonons, 

O = njf + njj. - 1 . (5) 

It should be noted that our definition of Xs( w ) differs 
from the usual one by a factor of four. We also calculate 
the quasiparticle weight z given by z = (1 — ReE'(O)) -1 
within the local approximation of the DMFT. 

The different regimes of the phase diagram of Fig. 1 
give rise to the following structure of the paper: In Sec. II 
we choose a fixed value of U and study dynamical prop- 
erties as a function of the phonon coupling g. We find 
continuous transitions from a metallic to a bipolaronic 
state for small values of U = and U = 1 (Sec. II A). 
For a larger value of U — 5, as discussed in Sec. II B, 
the transition to the bipolaronic state becomes discon- 
tinuous. Complcmcntarily, we look at the [/-dependence 
of dynamical properties in Sec. Ill keeping g fixed to 
g = 0.45. For this value of g the system can be in all 
three phases, depending on U . In Sec IV we study the 
properties of the model on the two polaronic lines as given 
by U e ft — and {n^n{) = 1/4 (dashed and dot-dashed 
lines in Fig. 1). We also show spectra for U e g = for 
different values of u>q. In Sec. V we give an overview of 
our results, and we present our conclusions in Sec. VI. 

II. DEPENDENCE ON g - PHONON DRIVEN 
PHASE TRANSITION 

In this section, we fix the value of the Hubbard inter- 
action U and study dynamical properties for a varying 
electron-phonon coupling g. We distinguish between the 
cases of weak and strong U. 

A. Weak U 

First of all, we look at the one-electron spectra for 
a fixed U = 1 as shown in Fig. 2. There is a strong 
similarity of these results with those of the pure Holstein 
model (U = 0) calculated previously 16 . At weak coupling 
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FIG. 2: One-particle spectral function for (7 = 1 and various 
values of g. The inset shows the variation of the quasiparticle 
weight z for U = 1 (solid line) and U = (dashed line). 



a narrow peak emerges at the Fermi energy. The top 
of this feature is rather flat and the imaginary part of 
the self energy has a shallow w 2 -dependence. This one 
would expect from lowest order perturbation theory in g, 
where the imaginary part of the self energy vanishes in 
the range \w\ < uj q for U — (see Eq. (19) in Ref. 24). 
Upon increasing g further, there is a rapid narrowing of 
this feature until it disappears at a critical coupling of 
g c w 0.47 and a gap opens. Similar to the pure Holstein 
model, there is no preformed gap in this case in contrast 
to the Mott transition in the Hubbard model. 

The metallic regime corresponds to a renormalized 
Fermi liquid. The quasiparticle weight z is shown as a 
function of g in the inset together with that for the pure 
Holstein model (U = 0). For U = 1, the quasiparticles 
are already slightly renormalized at g = 0. However, ini- 
tially the renormalization with g is weaker than in the 
case U = but more rapid on the approach to g 22 . 
As with the pure Holstein model, there is no evidence 
of multiple solutions near the phase transition. This is 
in contrast to the situation near the Mott transition for 
the Hubbard model, where there is a finite coexistence 
region. 

Figure 3 shows the spectra of the corresponding 
charge- and spin susceptibilities at U = 1 and {7 = 
for comparison. In the free system, x c (cj) = x s (c<j) as 
defined in Eq. (4), which can be seen in the lowest curves 
of the two upper panels. Looking at the case U = first 
(upper panel), we see the emergence and buildup of a 
low-energy peak in Xc and the suppression of weight in 
Xs for uj < W/2. In the bipolaronic state (g — 0.42), 
both charge and spin susceptibilities have the same gap, 
roughly twice that seen in the electronic spectrum. They 
have an identical peak above this gap. These can be seen 
in the dotted curves of Fig. 3, which have been scaled 
up by a factor of 10 to make these features visible. The 
coincidence of these peaks can be explained by the fact 
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FIG. 3: Low energy behaviour of the spectra of the charge 
(left) and spin-susceptibility (right) for the pure Holstein 
model (U = 0, top) and U = 1 (bottom). The values of g 
are g = 0.0 (bottom), g = 0.1, 0.2, 0.3, 0.37, 0.40, 0.42 (top) in 
the upper panel and g — 0.0 to 0.5 in steps of 0.1 in the lower 
panel. The dotted lines are scaled up by a factor of 10. 

that spin fluctuations necessarily require the breakup of 
a bipolaronic pair. Their weight is very small due to the 
strong bipolaronic binding. There is an additional very 
small peak in the charge susceptibility in the gap close 
to ujq from the residual couplings to the phonon mode. 
In the absence of charge ordering, there must be charge 
fluctuations at uj = which, however, cannot be observed 
numerically. 

For U = 1 (lower panel) we observe the same overall 
trends. However, starting from a reduced charge suscep- 
tibility due to the finite U, the peak in Xc develops more 
slowly. The spin susceptibility, which is enhanced due 
to the finite U, changes very little in the metallic phase. 
The low-energy spin fluctuations are almost completely 
suppressed with the emergence of the gap in the bipola- 
ronic state. 

The charge fluctuations can be directly related to the 
phonon propagator via the equation (see appendix A) 

d(uj) = do(w) + g 2 d (uj) Xc(u) <k>{u) , (6) 

where do(u>) — 1/(cj — ujq) is the free phonon propagator. 
However, use of this equation with an approximate form 
of Xc{uj) can lead to severe numerical errors, as discussed 
in appendix A. In Fig. 4, we give the results for the 
phonon spectra for U = and U = 1 calculated via the 
self energy 28 , similar to the procedure we use for the elec- 
tronic spectra 27 . The results for U — differ significantly 
from those published earlier 16 based on Eq. 6. 

In contrast to the earlier results, the main phonon peak 
is found initially to soften significantly in a similar way to 
the Migdal Elisahberg result. The two-peak structure de- 
velops only upon approaching the transition to the bipo- 
laronic state. In the gapped state, the high-frequency 
peak narrows and tends towards u>q. The same trend 
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FIG. 4: Spectral density of the phonon propagator for U = (left) and U = 1 (right) and various values of g. 



can be seen for U = 1 (right-hand plot). The slower ini- 
tital softening correlates with the suppression of charge 
fluctuations seen in x c ( w )- 

On the approach to the transition, significant nega- 
tive spectral weight for u> < can be seen in the spectra 
shown in Fig. 4. As a consequence of the spectral theo- 
rem, this relates directly to the average number of excited 
phonons in the ground state. The spectra of the Green's 
function D(uj) show similar features as d(u>), and inte- 
grated up to uj = yield the average lattice fluctuations 
as shown in Fig. 5. These results agree well with the di- 
rect evaluation of the expectation value ((£> + b^) 2 ) in the 
metallic phase. For low to intermediate values of g, the 
NRG results do not deviate significantly from the Mid- 
gal Eliashberg calculation and are significantly smaller 
than reported in Ref. 16 based on Eq. 6. On the ap- 



proach to the transition, there is a large increase in the 
lattice fluctuation. In the gapped phase, we observe a 
large difference between the direct calculation and the 
integration over the negative part of the phonon spectral 
density. One way of looking at it is that due to the lim- 
ited numerical resolution we miss the contribution from 
the peak at to = 0~ in the spectrum of D(u)). An other 
way of explaining this difference is that there is no dy- 
namics connecting the two degenerate ground states with 
(n) = and 2 in the gapped phase. These two ground 
states have nonvanishing average displacements ±xo, see 
appendix B, and integration over the phonon spectrum 
yields the fluctuations about ±x , i.e., ((x — x ) 2 ). The 
calculation of xq in one of the ground states quantita- 
tively explains the different values of the two methods in 
the gapped phase. 
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FIG. 5: Expectation value of the lattice fluctuations as a func- 
tion of the electron-phonon coupling g in the pure Holstein 
model (U = 0). In the metallic range (g < 0.4), the direct cal- 
culation (+) agrees well with the integration over the phonon 
Green's function (o) and charge susceptibility (x)). For de- 
tails see appendices A and B. The dot-dashed line shows the 
result from the Migdal Eliashberg calculation. 



B. Strong U 

Next we look at the transition to the bipolaronic state 
for a larger fixed value of U = 5. Figure 6 shows the 
corresponding one-particle electron spectral functions for 
various values of g. Their ^-dependence is in sharp con- 
trast to the small U case shown in Fig. 2. Initially for 
g = 0, we have the three-peak structure of a strongly 
correlated Hubbard model. With increasing g, the cen- 
tral resonance broadens slightly as is reflected in the in- 
creasing quasiparticle weight z (see inset). This can be 
ascribed to the partial cancellation of the Hubbard repul- 
sion by phonon mediated electron-electron attraction, see 
Eq. (2). The effect, however, is much weaker than one 
would deduce from U e g. The position of the high-energy 
Hubbard bands is virtually unaffected by the phonons, 
as expected from Eq. (2) evaluated at u) « ±£7/2. When 
approaching the critical coupling g c rs 0.72, the quasi- 
particle peak disappears abruptly and z jumps to zero, 
indicating a discontinuous transition to the bipolaronic 
phase. The double occupancy, shown in the inset to 
Fig. 6, is much less than a quarter in the metallic state 
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FIG. 6: One-particle spectral function for (7 = 5 and various 
values of g. The central peak broadens with increasing g and 
at g ~ 0.71 disappears in a first order phase transition to the 
bipolaronic state. The inset shows the quasiparticle weight z 
(o) and the double occupancy d = {n<\n{) (o) as functions of 
9- 



FIG. 8: Spectral density of the phonon propagator for (7 = 5 
and various values of g. 

low-energy shoulder. After the discontinuous transition 
to the bipolaronic phase, the mode hardens back to loo 
and narrows. 



but makes a sudden jump to d w 1/2 on the point of the 
phase transition. A coexistence of metallic and gapped 
solutions is found in the narrow range 0.71 < g < 0.72. 

The spectra of the charge and spin susceptibilities are 
shown in Fig. 7. They display similar trends, but more 
pronounced, than the corresponding results for U = 1 
(lower panels of Fig. 3). The low-energy peak in Xc( w ) 
only becomes clearly visible when g reaches the value 
of g w 0.65. In the bipolaronic phase, the peaks above 
the gap in the charge and spin susceptibilities are not 
identical as in the U = case, but very similar. However, 
the gaps are too large to be visible on the w-scale of Fig. 7. 
The only peak visible, when enhanced by a factor of 20, 
is the peak in the charge susceptibility within the gap 
near uiq. 

In Fig. 8, the spectra of the corresponding phonon 
propagators are plotted. Over a large initial range of 
g < 0.6 there is very little softening, corresponding to a 
complete suppression of charge fluctuations. When ap- 
proaching g c , the phonon mode softens and develops a 
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FIG. 7: Low-energy behaviour of the spectra of the charge 
(left) and spin-susceptibility (right) for (7 = 5. The values of 
g correspond to those in Fig. 8. Note the different scales of 
the j/-axes. Dotted lines have been scaled up by a factor of 
10. 



III. DEPENDENCE ON U - MOTT 
TRANSITION 

In this section we will fix the electron-phonon coupling 
at the value g = 0.45 and study the dependence of dy- 
namic response functions on the Hubbard repulsion U. 
The value of g = 0.45 is chosen because all three phases 
can be found for it, depending on U (see Fig. 1). 

Figure 9 shows the one-electron spectral function in the 
vicinity of both phase transitions. The right-hand plot 
for large values of U is very similar to what is observed 
in the Mott transition of the pure Hubbard model. The 
lower and the upper Hubbard peak at ±(7/2 are well 
developed and move to higher energies as U increases. 
At the same time, the central peak narrows and vanishes 
at the Mott transition. As for the pure Hubbard model, 
there is a preformed gap in the one-electron spectrum 
and a broad region of numerical coexistance of metallic 
and insulating solutions. The electronic spectra are very 
similar to those of the pure Hubbard model near the Mott 
transition. The phonons do not seem to alter much the 
picture of the phase transition. 

The transition from the bipolaronic to the metallic 
state is completely different, as can be seen in the left- 
hand plot of Fig. 9. The phase transition occurs close 
to (7 = 1 and the spectra show a similar behaviour as 
discussed in Sec. II A, where g is varied. In the bipola- 
ronic state, with the initial increase of (7, the two polaron 
bands move towards the Fermi level. On the metallic 
side, just after the transition, we see a very sharp reso- 
nance which broadens when U is increased further. There 
is no signature of a preformed gap once the system is 
metallic. Again, there is no evidence of a significant co- 
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FIG. 9: One-particle spectral function for g = 0.45 and values of U close to the transition from the bipolaronic to the metallic 
state (left) and U close to the transition to the Mott insulator (right). 



existence region. 

The spectra of the charge and spin susceptibilities for 
g = 0.45 are shown in Fig. 10. We clearly see that the 
low-energy features of the two susceptibilities signal the 
appearance of two different instabilities of the system, 
depending on U. In the bipolaronic phase, both spectra 
are gapped and we see only the peak near u>o in the charge 
susceptibility. Upon entering the metallic phase, the zero 
frequency peak in the charge susceptibility moves to finite 
uj. This peak loses weight and slightly disperses to higher 
energies as the system becomes more metallic. In this 
small U regime, spin fluctuations are suppressed. 

The converse occurs when approaching the Mott tran- 
sition at U ~ 6.2. The charge fluctuations are suppressed 
and peak in the spin susceptibility builds up, signalling 
instability towards antiferromagnetic ordering. In the 
Mott insulating state, both susceptibilities are gapped. 
The spin fluctuations above the gap (beyond the scale of 
Fig. 10) must largely be due to the charge fluctuations, 
as in the pure Hubbard case. However these peaks are 
not identical, in contrast to the bipolaronic gapped state. 
The difference in the peaks can be interpreted as a result 
of the spin correlations induced by U . However, due to 
the phonon coupling the spin fluctuation peak is reduced 
compared to the pure Hubbard case. 

Figure 11 displays the phonon propagator for this scan. 
In the bipolaronic phase at U = 0, electrons and phonons 
are almost decoupled, and the phonon spectrum shows 
only one peak just below u>o — 0.2. When increasing U, 
this mode softens and splits into two peaks close to the 
phase transition. On the metallic side, the low-energy 
peak survives and continuously hardens back to ujq, as 
the system approaches the Mott transition. There is no 
signature of the Mott transition in the phonon spectrum. 
The effect of an increasing U is just to suppress contin- 
uously the charge fluctuations which results in a decou- 
pling of electrons and phonons which drives the harden- 
ing of the phonon peak. 



Another feature of the suppression of charge fluctua- 
tions is the narrowing of the phonon peak with increase 
of Z7, as also observed in Ref. 18. In our case, however, 
the shift of the peak appears to be more marked. The 
narrowing can also be seen in comparing the curves with 
the same (low) value of g in Figs. 4 and 8. 

Further insight into the nature of these transitions can 
be gained by looking at the behaviour of the higher order 
electron Green's function 

F„(u) = ((cAc* (7) 

This Green's function is needed for the calculation of 
the self energy in the NRG procedure 27 . The integrated 
spectral weight for — oo < uj < equals the double oc- 
cupancy (n^ni) 27 . Moreover, the total spectral weight 
of F„(w) yields the average density (n ff ), which should 
be (rtf) = (rij) = 0.5 in the particle- hole symmetric case 
treated here. 

Figure 12 shows the spectral function of F(u>) for 
g = 0.45. In the bipolaronic phase, almost all of the 
weight is located below the Fermi energy signalling a high 
value of the double occupancy (see inset). As we enter 
the metallic phase, a large part of this weight is trans- 
ferred to higher energies. In the Mott insulator, virtually 
no weight is left at u> < and the double occupancy is 
completely suppressed. 

The inset to Fig. 12 also shows the total spectral weight 
of F(u>), which should evaluate to 1/2. We see that this is 
indeed more or less the case, except for values of U close 
to the phase transitions U m 1 and U » 6. For these 
values, the truncation of the Hilbert space in the NRG 
introduces the most significant errors. There we also find 
the strongest discrepancies between the direct calculation 
of the double occupancy from the ground state expecta- 
tion value (marked with x ) and the indirect, less accurate 
calculation via F(u>) (marked with o). 
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FIG. 10: Spectra of the charge (left) and spin-susceptibility (right) for g — 0.45 and various values of U (same as in Fig. 11). 
The charge susceptibility reflects the transition from a bipolaronic state to a metal, whereas the spin susceptibility signals the 
transition to the Mott insulator. The dotted lines are scaled up by a factor of 100 to show the peak at ui = t^o in the bipolaronic 
phase. 
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FIG. 11: Spectral density of the phonon propagator for 
g — 0.45 and various values of U. The transition from the 
bipolaronic to the metallic state is clearly visible by the two- 
peak structure whereas the transition to the Mott insulator 
does not obviously affect the phonon propagator. 



IV. POLARONIC LINES 




FIG. 12: Spectral density of the higher electronic Green's 
function F{uj). The weight for u> < indicates the degree of 
double occupancy; g — 0.45 and various values of U. The inset 
shows the double occupancy calculated directly (x) and via 
the spectral theorem from F(u) (o). Circles indicate the total 
spectral weight of F(lj) which should evaluate to n — 1/2 at 
half filling. The discrepancies are most pronounced near the 
phase transitions at U as 1 and U ~ 6. 



In this section we discuss the dynamical properties 
along two particular lines in the g—U plane (see Fig. 1). 
The first of these lines is given by U c s = U + 2g 2 /ujq = 
where one might naively expect free quasiparticles close 
to the Fermi surface. The other line is the location of 
points where the double occupancy (n^n^) = 1/4, as in 
the free system. 

Only for very large values of luq can the Hubbard term 
be completely cancelled by the phonon-mediated attrac- 
tive term (see Eq. 2), in which case both polaronic lines 
would coincide. We investigate these lines to understand 
the degree of compensation of these competing interac- 



tions for small phonon frequency ujq = 0.2. 



A. Polaronic line {7 e ff = 

Figure 13 shows the one-electron spectral functions 
for various values of g and U such that U e g = 0. For 
small values of U (and g), we see the sharp feature in 
the electronic spectrum as discussed in Sec. II. How- 
ever, upon increasing U further, Hubbard bands develop 
and the central feature becomes the central quasiparti- 
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FIG. 13: One-particle spectral function on the line U e s — 
at u>o = 0.2 (dashed line in Fig. 1). The inset shows the 
variation of the quasiparticle weight z and double occu- 
pancy d = (n T n|>. 
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FIG. 14: Spectral density of the phonon propagator on various 
points on the polaronic line U g = 0. The inset shows the 
expectation value of the lattice fluctuations as a function of 
U. 



cle peak similar to the pure Hubbard model for large 
U(< U c ). This indicates that phonons play their as- 
sumed role of compensating the Hubbard repulsion only 
to a very limited extent. The only remarkable difference 
to the pure Hubbard model is that on the line U e s = 
the phase transition is shifted towards a higher value of 
U c ~ 6.5, as compared to U c — 5.88 for g — 0. The inset 
to Fig. 13 shows that the quasiparticle weight decreases 
steadily from z = 1 to zero and the double occupancy 
from (n^rii) — 1/4 to a small but finite value. 

The phonon spectra, as displayed in Fig. 14, are always 
made up of only one single peak that initially moves to- 
wards lower energies as more phonons are excited and 
lattice fluctuations increase (see inset). At U ~ 3 this 
trend reverses and the phonon peak hardens back when 
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FIG. 15: One-particle spectral function on the line (n^ny) — 
1/4. (dot-dashed line in Fig. 1). The inset shows the variation 
of the quasiparticle weight z along this line. 



approaching the phase transition to the gapped state. Si- 
multaneously, the number of excited phonons decreases 
to zero and the lattice fluctuations reach the value of the 
quantum mechanical zero point fluctuations. The reason 
for this non-monotonic behaviour of the lattice fluctua- 
tions as function of g is not obvious. The maximum in 
the lattice fluctuations as shown in the inset correlates 
with the value of U where the transition to the bipola- 
ronic state changes from second order (for U < 3) to first 
order (for U > 3). This could be a possible explanation 
as one would expect to see stronger lattice fluctuations 
as one aproaches a second-order transition to the bipola- 
ronic state. 



B. Locus of points with (n-\ny) = 1/4 

The line in the phase diagram where (n-j-n^) = 0.25 
starts at the uncorrelated system (g = U = 0) and ends 
at 17 w 3 where it merges with the phase boundary to 
the bipolaronic state. As a consequence, for larger values 
of U, the double occupancy has a jump at the transition 
from the metallic state ((n|n^) < 0.25) to the bipolaronic 
state ((n^ni) fa 0.5) (see Fig. 3 in Ref. 22 and inset of 
Fig. 6). Therefore this transition has to be first order in 
this range. 

Figure 15 shows the one-electron spectral functions 
along this line. We observe the appearance of the nar- 
row feature at the Fermi energy that is a typical for 
the electron-phonon coupling. The band develops broad 
shoulders, but they never separate into distinct sub- 
bands. 

The phonon propagator, as displayed in Fig. 16 shows a 
remarkable softening of the original peak and the build- 
up of negative spectral weight. As a consequence, the 
lattice fluctuations increase almost linearly with U and 
are quite pronounced for U — > 3. 
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FIG. 16: Spectral density of the phonon propagator on var- 
ious points on the line (n^ni) = 1/4. The inset shows the 
expectation value of the lattice fluctuations as a function of 
U along this line. 



C. Dependence on luq 

So far we have worked with a fixed value of uj = 0.2. 
However, the degree of compensation of the competing 
interactions is dependent on ujq. Here we return to the 
case U c s — and study the variation of loq with fixed 
U = 5 and g = ^/Uuq/2. 

In Fig. 17 we plot the one-electron spectral density for 
a range of values of u>q. For the largest value of u>o = 16, 
the spectrum is virtually the same as that of the free 
system, because we have almost complete compensation. 
We now examine the behaviour as we progressively re- 
duce u>q. First the high-energy phonon subbands gain 
weight and move towards the Fermi energy. Their po- 
sition is roughly given by u)q. There is a commensurate 
reduction of the width of the central peak. For u!q = 2 the 
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FIG. 17: One-particle spectral function for U e s = 0,U — 5 
and various values of u>o ■ 



FIG. 18: Spectral density of the phonon propagator for [/ e ff 
and [7 = 5 for various values of ojq ■ 



phonon subbands are visible as shoulders of the central 
peak. Upon further reduction of u>q, these shoulders move 
back towards higher energies again and, for uj < 0.5, be- 
come the Hubbard bands located at uj ss ±U/2 = ±2.5. 
For the smallest value of luq = 0.2, the spectrum is iden- 
tical to the one discussed in Sec. IV A. 

In Fig. 18 we plot the spectra of the phonon propagator 
as a function of the relative energy scale lu/loq. For large 
luq = 16, where the compensation is almost complete, we 
see the narrow phonon peak at loq. In addition we see a 
small feature near u = 0. As ujq is decreased, the low- 
energy feature gains weight and the mode at luq broadens 
and moves to slightly higher energies, as one would ex- 
pect when too approaches but lies above the upper band 
edge. At ojq = 4, the two features have approximately 
equal weight. For smaller values of luq the upper peak 
has almost zero weight. The remaining peak is the soft 
mode discussed in the previous sections. 

In Fig. 19 we illustrate how different static quantities 
reflect the trend of the Holstein-Hubbard model towards 
a free system in the limit of large u>o- We plot z and 
Ain^ni) as a function of u>q. The double occupation 
value converges relatively rapidly to that of the free state. 
The quasiparticle weight, however, converges rather more 
slowly. In fact, only for ujq 3> W — 4 is it approaching 
the uncorrelated limit z = 1. 

In the same Fig. 19 we also plot the lattice fluctuations 
(x 2 ) = {(b + b^) 2 ) / (2u>o)- This curve goes monotonically 
to zero as ojq — > oo asymptotically behaving as 1/loq for 
large ojo- However, if we plot the rescaled values ((b + 
b^) 2 }, we find a a maximum which coincides with the 
region of maximum softening of the phonon mode. For 
larger values of u>o this quantity decreases to the quantum 
limit of (0 + 6t) 2 ) = 1. 

One would expect the phonons to decouple from the 
electrons in the limit of large luq- It is interesting to ob- 
serve, however, that 1+4(6^6} (dot-dashed line in Fig. 19) 
fits well the high-^o behaviour of the lattice fluctuations 
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FIG. 19: Variation of the quasiparticle weight z (o), the 
rescaled double occupancy 4/n-fnj.) (o), and lattice fluctua- 
tions (a; 2 ) (x) as a function of wo for U e s = 0. The dot-dashed 
line represents the number of excited phonons 1 + i(b'b). 

which indicates that the ground state of the phonons in 
this limit is a coherent state rather than an eigenstate of 
the free phonon Hamiltonian Lo^wb . 

V. OVERVIEW 

For an overview we will find it useful to classify the 
metallic state in the phase diagram into three regions, 
separated by the two polaronic lines, In each of these we 
find qualitatively different behaviour. 

We consider first of all the region below the line U c s = 
0. Here U c s > 0. One interesting question is how the cou- 
pling to the phonons affects the transition to the Mott 
insulating state. Our results indicate that this transition 
is always very similar to that found in the pure Hubbard 
model: the appearance of a preformed gap, a rather large 
coexistence region between insulating and metallic solu- 
tions and the quasiparticle weight continuously going to 
zero. There is only a small increase in the critical value 
of U which is found empirically to depend upon g to good 
approximation as 

U c « C/ c , Hu bb + 0.8 x g 2 . 

The small influence of the phonon coupling can be read- 
ily explained by suppression of charge fluctuations for 
large U. As the electron-phonon coupling in the Holstein- 
Hubbard model is to the charge fluctuations, the suppres- 
sion of these renders the coupling ineffective. In fact, the 
whole region U c s > appears to be dominated by the 
Hubbard term. For w > wo the second term in the re- 
tarded interaction as defined in Eq. 2 is negligible and 
the effective interaction is essentially given by U. Also 
on the lowest energy scales u> <C oJq, U e g(uj) > 0. 

Next we look at the complementary region defined by 
(n-[Tii) > 1/4 which is dominated by the coupling to 



phonons. Here the similarity is with the pure Holstein 
model. The metal to bipolaronic transition takes place 
without a preformed gap and with the characteristic soft- 
ening of the phonon mode. The transition in this region 
appears also to be continuous with no significant coex- 
istence region. The main effect of the Hubbard U is to 
push the phase transition and its precursors to somewhat 
larger values of g. 

Finally, in the region bounded by the two curves 
U c s — and (n^n^) — 1/4 there is more complex in- 
terplay of the two interactions, as can be seen in our re- 
sults on the response functions. Beyond the point where 
the line (n-|-n^) = 1/4 merges into the phase boundary 
the transition to the bipolaronic state is first order and 
thus qualitatively different from that discussed in the last 
paragraph. Along the lower boundary line of this region, 
i.e., for U e g = 0, the effective interaction is still repul- 
sive for w > wo. The question arises whether on the 
lowest energy scale there is a complete cancellation of 
the two competing interactions. This question can be 
investigated by examining the lowest-lying one- and two- 
particle excitations from the interacting ground state. 
The calculation of the effective interaction U e s between 
the renormalized quasiparticles has been considered in 
detail for impurity models in Ref. 29. A similar analysis 
is possible using the effective impurity model within the 
DMFT. From Ref. 29, 

C/ eff oc lim (E pp (N)-2E p (N))A ( - N ~ 1 V 2 , (8) 

where E P (N) (E PP (N)) is the energy of the lowest- lying 
single-particle (particle-particle) excitation in the N-th 
iteration of the NRG procedure and A > 1 is the NRG 
discretization parameter. As for the impurity model, we 
find that, within the numerical accuracy of the NRG, U e g 
vanishes on the line U e g = and is positive (negative) 
below (above) this line. Therefore, along this line we have 
a system of renormalized (z < 1) but non-interacting 
quasiparticles. On the high-energy scale lo > luq, the 
system is still dominated by U. This can be seen by the 
formation of the Hubbard bands as shown in Fig. 13 and 
the suppression of the double occupancy. 

VI. CONCLUSIONS 

In this study we have investigated the particle-hole 
symmetric Holstein-Hubbard model. We have used the 
dynamical mean field theory in combination with the 
numerical renormalization group to calculate dynami- 
cal correlation functions for the full Hamiltonian with 
quantum phonons. This non-perturbative approach al- 
lows us to access all parameter regimes of the model 
at zero temperature. A non-perturbative technique is 
indeed necessary because perturbative methods such as 
Midgal-Eliashberg approach for the Holstein model are 
known to break down 15 in the strong-coupling limit. 

We find three regions with qualitatively different be- 
haviour. In the region of the Mott metal insulator transi- 
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tion we find the phonon effects are largely suppressed by 
the on-site repulsion U. This conclusion may be specific 
for the non-degenerate Holstcin-Hubbard model where 
the coupling is purely to the local charge density. A 
system with orbital degeneracy and Jahn- Teller phonons 
shows 20 strong phonon effects even for strong U. In our 
case of the non-degenerate Holstein-Hubbard model, the 
electronic effects on phonons yield only a modest soft- 
ening of the phonon mode in the metallic regime. This 
softening disappears in the Mott insulator. 

Weak electron-electron interactions seem to have little 
effect other than to delay the onset of the transition to the 
bipolaronic phase. For larger values of U the transition 
changes from continuous to discontinuous. For small val- 
ues of U we see a complete softening of the phonon peak 
when approaching the transition. The softening is much 
less pronounced for larger U. 

The line U e g = in the metallic phase does seem to 
have some significance in that for J7 off > the spin sus- 
ceptibility has a dominant low-energy peak. For U e s < 
the dominant low-energy peak is in the charge suscepti- 
bility. The line U c s = also appears to be significant 
for the low-energy Fermi liquid behaviour as the on-site 
quasiparticle interaction changes sign on or close to this 
line. This is not obvious as one might have expected 
that on the low energy scale of the electron-phonon in- 
teraction, the electron-electron interaction is no longer 
given by the bare U but some renormalizcd value U . This 
would imply that one should use U in Eq. 2 on the energy 
scale w <C u)q, so that the change of sign of the quasipar- 
ticle interaction U c s, as defined in Eq. (8), would occur 
when U — 2g 2 / ujq = 0. However, our results indicate that 
this is not the case and that U e g vanishes on or close 
to the line U e g = U — 2g /uo = 0. A similar situation 
was found in the case the Anderson-Holstcin impurity 
model 29 . 

The Holstcin-Hubbard model is very rich and shows 
also diverse forms of behaviour that could not be ad- 
dressed in this paper. Apart from the results pre- 
sented here, this model is known to exhibit various types 
of symmetry-breaking phases 17 . We intend to extend 
our calculations to investigate the competition between 
charge order, antiferromagnetism, superconductivity and 
the effects of doping 30 . This should enable us to make 
direct contact with experimental results on such com- 
pounds as the fullerides, where this model with the as- 
sumption of the coupling of the electron density to the 
local phonon mode should be directly applicable. For 
applications to Jahn- Teller systems, such as in the man- 
ganites, the model would have to be extended to include 
degenerate orbitals. 
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APPENDIX A: PHONON PROPAGATOR AND 

XcK) = o 

Here we reconsider the relation between the phonon 
propagator and the charge susceptibility (Eq. (6)) for the 
Holstein-Hubbard model and comment on pitfalls in its 
exploitation. 

Equations of motion: For the derivation of Eq. (6) , we 
can use the standard equation of motion for the Fourier 
transform of the double-time Green's function, 

cj ((A ; B)) u = ([A, B] v ) + (([A, H}_ ; B)) ^ (Al) 

For the boson Green's function d(w) we take A — b i , B = 
b\, and r\ = — 1. For D(lo) we take A = B = b i + b\ and 
T] = — 1. The following procedure works for both d(u>) 
and D(lj), but we will use the former in what follows. 
The equation of motion gives 

{w-u )d(u>) = l + g((d i ;b\)) u , (A2) 

where Oj = n il7 — 1. Taking the equation of motion 
for the right-hand operator, we have 



(A3) 



Hence the result (do(w) = (oj — ujq) 1 is the free phonon 
propagator), 

d(Lu) = do(w) + g 2 d (Lu) {(Oi ; 6i)) u do(f>) ( A4 ) 

which is Eq. (6) linking the phonon propagator d(u>) with 
the charge susceptibility Xc(u) = ((Oi ; O,)) . 

As a Green's function, the phonon propagator has a 
series of poles of first order, but none of second order. 
For this to be the case, Eq. (A4) tells us that the charge 
susceptibility has to vanish as 



Xc{u) ~ (lj - w ) for lo — > loq 



(A5) 



in order to ensure that the pole of d(u>) at lu = ujq is first 
order. This is also evident if we express x c (w) in terms 
of the irreducible particle-hole bubble U(u). It takes the 
form 



Xc(w) = 



1 - gm(u)D (cu) 



(A6) 



with the non- interacting phonon propagator D (uj) = 
2u) /(u> 2 — u)q). Since D (tu) diverges at w = uj , we 
have Xc(^o) = 0. 

Spectral Density: We will derive an expression for the 
spectral density of the phonon propagator, i.e., for 



PdioS) = lim d(u> + iS) 



(A7) 



12 



that implicitly includes the fact that Xc(wo) = 0. 
Consider the system with a finite basis such that, 



Xc 



M = E 



a, 



a; — w, 



(A8) 



for u> > 0, where {w^} are the discrete set of excitations 
in the charge density response function. Then 



(lu — luq + i5) 



Xc(u + i5) = -Xc(u + id)- 



1 



dto (to — ujq + i5) 
(A9) 

The contributions to the spectral density Pd{u) arise 
solely from the poles at ui = loq and to = uij. These 
do not coincide because Xc(wo) = 0- The contribution 
to Pd(u) from the poles at to — uij are straightforward to 
evaluate and give 

The contribution from the pole at lu = loq is a little trick- 
ier to evaluate. The term in g 2 contributes 



-9 2 Xc(.u)5'(u - w ) 



(All) 



where 5'(lo — loq) denotes the derivative of the delta func- 
tion. However, we have 

Xc(u)S'(uj - w ) = -x' c (uo)o~(u - w ) + Xc(wo)^'(w - oj ) 

(A12) 

where x' c (w) is the derivative of Xc(w) and is given by 



met. Therefore Xc(w) has no root exactly at to = luq. 
This leads to a double pole of cL(lo) at ui = loq and to 
a distorted result for the phonon propagator at energies 
close to loq , when derived from Eq. (A4). Lattice fluctu- 
ation, as calculated by integrating D(uj) from Eq. (A4) 
are strongly overestimated for the same reason. 



APPENDIX B: FORMULA FOR AVERAGE 
DISPLACEMENT (x): 

In this appendix we derive an expression linking the 
average local displacement (x^ ~ (b\ + b t ) to the lo- 
cal electron density. For this derivation, we modify the 
Hamiltonian by coupling an extra term linear in Xi at 
each site. The extra term in the Hamiltonian reads 



(Bl) 



with the coupling a = {en}. Upon differentiating the 
expression for the partition function Z with respect to 
a;, we obtain 



{b * +b * ) = -z]3 0^ 



a=0 



We then apply a canonical transformation to the Hamil- 
tonian, H = U~ 1 HU ', with U given by 

(] = TT e -(« i +s6 i )(6l-6 i )/a;o _ 



X, 



(W ~ U)j) 2 



(A13) 



Using Xc(wo) = and collecting all the terms together 
we obtain 



Pd 



M= i-.9 2 E 



f ( w o-Wj) 2 



S(lo - ljq) 



(A14) 



+ q 2 V S(lo - to., 

9 ^(uq-uJj) 2 y J 



This expression can be used to calculate the spectral 
density of the phonon propagator Pd(u) from the peaks 
{aj,(jj} of the charge susceptibility Xc(w). In fact, the 
numerical evaluation of Eq. (A14) yields the correct re- 
sult for the phonon propagator and the lattice fluctua- 
tions, as shown in Fig. 5 for the case of the pure Holstcin 
model. The condition Xc(wo) = has, of course, been 
used explicitly in the derivation of Eq. (A14). 

The problem of the direct evaluation of Eq. (A4) for 
the phonon propagator is the numerical error in Xc(w). 
In fact, the truncation of the Hilbert space in the NRG 
procedure entails that the condition (A5) is not exactly 



This is a displaced oscillator transformation at each lat- 
tice site. The phonon and electron operators are trans- 
formed as 



b l = U- 1 b l U = b l -^-^-6 i , 
wo w 

After some algebra we find that the terms depending on 
a in the transformed Hamiltonian H = U~ l HU read 



w w / 



(B2) 



The partition function Z is, of course, not changed by 
the canonical transformation of the Hamiltonian. Upon 
expressing Z in terms of H , differentiating with respect 
to on and then putting a = we find the result 



1 



which relates the average local displacement to the local 
average electron density. 
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